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The Global Positioning System (GPS) has become a standard 
method for low cost onboard satellite orbit determination. The 
use of a GPS receiver as an attitude and rate sensor has also 
been developed in the recent past. Additionally, focus has 
been given to attitude and orbit estimation using the 
magnetometer, a low cost, reliable sensor. Combining 
measurements from both GPS and a magnetometer can provide 
a robust navigation system that takes advantage of the 
estimation qualities of both measurements. Ultimately, a low 
cost, accurate navigation system can result, potentially 
eliminating the need for more costly sensors, including 
gyroscopes. This work presents the development of a 
technique to eliminate numerical differentiation of the GPS 
phase measurements and also compares the use of one versus 
two GPS satellites. 

INTRODUCTION 

This work presents the testing of a unified navigation algorithm that produces 
estimates of attitude, angular rate, position, and velocity for a low earth orbit (LEO) 
spacecraft. The system relies on GPS phase and range data as well as magnetometer 
data. The algorithm used is an extended Kalman filter (EKF), blended with a 
‘pseudo-linear’ Kalman filter algorithm, developed to provide LEO attitude, orbit, 
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and rate estimates using magnetometer and sun sensor data 1 . For many LEO 
spacecraft the sun data is available during only a portion of the orbit. However, GPS 
data is available continuously throughout the orbit. GPS can produce accurate orbit 
estimates and combining GPS and magnetometer data improves the attitude and rate 
estimates 2 . The magnetometer based EKF can converge from large initial errors in 
position, velocity, and attitude 3 . Combining the magnetometer and GPS data into a 
single EKF provides a more robust and accurate system. 

The algorithm is based on the EKF of Ref. 1. The GPS measurement models for 
phase and range are incorporated into the existing structure of the filter. The original 
EKF produced the orbit estimates in terms of Keplerian elements. Due to the nature 
of the GPS measurements and ease of computation, the orbit estimates are converted 
to the Cartesian position and velocity. The measurement model for the 
magnetometer is adjusted for this change in the state. 

The estimation of the angular velocity, or rate, typically involves numerical 
differentiation that introduces noise into the update of the rate estimate. An 
alternative approach, called the ‘estimation approach’, avoids differentiation by 
incorporating the kinematics into the filter 4 . This method, applied to the 

magnetometer, is reported on in Ref. 5. 

In the present work, the ‘estimation approach’ is also applied to the GPS phase 
measurements. Further testing with noisy sensor data is conducted on the EKF. 
Comparisons are made between the ‘estimation approach’ and the ‘derivative 
approach’ (requires numerical differentiation in the rate estimation ). Additionally, 
comparisons are made between estimates computed with a single GPS update and 
two GPS updates to determine if the additional complexity of two GPS updates 
provides a significant increase in the accuracy of the estimates. 

ALGORITHM DEVELOPMENT 

The GPS/magnetometer navigation system is based on an extended Kalman filter 
algorithm, blended with a ‘pseudo-linear’ Kalman filter originally presented in Ref. 
1. (Further developments are presented in Ref. 5.) In the present work, the 
development of the ‘estimation approach’ as applied to the GPS phase measurements 
is presented. The EKF algorithm is summarized only. The general EKF equations 
are: 


X(t) = f(X(t),t) + w(t) 

(i) 

k+i = ii k+ i (X(^k+i )) + iik+i 

(2) 


The general equations for the ‘pseudo-linear’ Kalman filter are: 

X(t) = F(X(t k ))X(t) + w(t) 

2k+i : H k+1 X k+l + n k+l 


The output filter state is given as 


( 3 ) 

( 4 ) 


( 5 ) 


x T = [E T ,Y T ,c T , a T ,oi T ,b T ,z T ] 

Where R and V are the spacecraft position and velocity vectors, respectively, C is a 
vector of GPS receiver clock errors, q is the attitude quaternion, © is the rotation rate, 
b is the magnetic field in body coordinates (a result of using the estimation 
approach’ with the magnetometer data 5 ), z is the GPS phase measurement which 
contains either the phase measurement states for one GPS observation or for two. 
The state, z, is developed in this work. 

The filter processes data from the magnetometer, a 2-baseline GPS receiver, and 
reaction wheels (used in Euler’s equations which describe the rotational dynamics of 
a spacecraft). As mentioned, the GPS data consists of range and phase 
measurements. Nominally, the filter first propagates the state estimate and filter 
covariance and then performs an update with the magnetometer data followed by an 
update with the GPS data. Data from either a single GPS satellite or from two GPS 
satellites are included in each update. The propagation of the state estimate and the 
filter covariance will be summarized first. The description of the state and 
covariance update will follow. 

Filter Propagation 


The dynamics for the position and velocity, the quaternion, and co, b, and z are 
treated as uncoupled. The position and velocity are propagated numerically with a 
two body gravity model which includes the j2 term. The clock errors are propagated 
according to a two state random-process model 6 . The quaternion is propagated 
numerically using the kinematics equation 7 . 

The rotation rate is propagated using Euler’s equation. The states b and z are 
propagated with the equation expressing the time rate of change of a vector with 
respect to a rotating coordinate system. The equations for the rate and b are 
presented in Ref. 5. Following is the development of the equation expressing the 
time rate of change of the elements of z. 

The expression for a single GPS phase measurement is given as 

Zij = • Sj (6) 


where: z,j = phase measurement 

5; = unit vector in direction of baseline i (i=l,2) 

Sj= unit vector in direction of a GPS satellite (j=l or 2) 

Equation (6) can also be written as (dropping the subscripts) 

z = a T s b = a T Ds, (7) 

where: a = baseline direction expressed in body coordinates 

D = direction cosine matrix, transforms from inertial to body coordinates 


si = satellite to GPS direction expressed in inertial coordinates 

Taking the derivative of (7) results in the following 

z = a T Ds, +a T Ds, (8) 

The time derivative of the direction cosine matrix can be written as 

D = -[cox]D (9) 

where: |rax] = cross product matrix of the rotation rate 

Substituting (9) into (8) and reordering the vectors gives 

z = a T [s b x]o> + a T D|, (10) 

where: [s b x] = [ds, x] = cross product matrix of the GPS direction unit vector in 
body coordinates 

D = current estimate of the direction cosine matrix which transforms from 
inertial to body coordinates 

s, = derivative of the GPS direction unit vector in inertial coordinates, this is 
not computed numerically but is developed from the relative position and 
velocity vectors (see Appendix) 

The vector z is formed from the individual phase measurements, z,j. (either two 
components for the two baselines and one GPS satellite or four components for the 
two baselines and two GPS satellites). Equation (10) is then written in the general 
form as 


z = F z co + u 


( 11 ) 


Equation (1 1) is used to propagate z numerically. 

In general, the linearized error dynamics are given as 

x(t) = Fx(t) (12) 

where: F= d(/dX 

The full dynamics matrix, F, is given in References 1 and 5 for all the states except 
the phase measurement states, z. The matrix F z from (11) is augmented into the 
matrix F. 

The filter covariance matrix is propagated according to 

p k+ , (-) = O k (X k (+))P k (+)O k T (X k (+)) + Q k 


(13) 


where: T> = the state transition matrix = e FAI 

Q = matrix of standard deviations of the white noise processes driving each of 
the states 

The filter uses the ‘multiplicative’ approach in the attitude estimation 3 , therefore the 
error state for the attitude is a vector of attitude errors. 

Filter Update 

The EKF state and covariance update equations are given as 

X k+1 (+) = X k+1 (-) + K k+1 [y k+1 - h k+ , (X k+1 (-)] (14) 

P k+1 (+) = [I - K k+1 H k+1 (X k+1 (-)]P k+1 (-)[I - K k+1 H k+1 (X k+1 (-)] T + K k+l RK k+l (15) 

The Kalman gain is computed according to 

K k+1 = Pk+i(~)H k+1 T (X k (-))[H k+1 (X k (-))P k+1 (-)H k+| T (X k (-)) + R k+1 ] _l (16) 


where H = 5h / 9X 

The measurement matrix, H, is partitioned as H = [H or bjt, H a ttitude> H ra te,b,z]- The 
individual measurement matrices are defined according to the two measurements 
used, namely the magnetometer measurement and the GPS phase and range 
measurements. The development of the magnetometer and GPS measurement 
matrices are given in Ref. 5, only the development pertinent to the additional states, 
z, is included here. The development of the GPS measurement matrix for the 
position and velocity states, presented in Ref. 5, includes both range and range-rate 
type measurements. In the current version of the EKF, only range measurements are 
included. 

GPS Phase Measurement State Update 

According to the development of the ‘estimation approach’ given in Ref. 4, the 
measured phase z m , is related trivially to the rate and the state z as 


z m =[o o i; 


CO 

b 

z 


+ = H rate b.,X,ate.b.z + I! 


where: p = added white noise 


(17) 


The innovations are then computed as 


and are used in the update equations for a linear Kalman filter. The H ra te,b,z matrix 
replaces the portion of the GPS measurement matrix developed for the ‘derivative 
approach’ given in Ref. 5. 

The use of the ‘estimation approach’ for GPS requires the initialization of the z 
states whenever a new GPS satellite is being observed. The states are initialized a 
priori with the first phase observations. The covariance is initialized with the a priori 
values and the correlation terms are set to zero. 


SIMULATION AND TEST CONDITIONS 

The test scenario consisted of an orbital ephemeris of the Rossi X-Ray Timing 
Explorer (RXTE) spacecraft and simulated orbits for the GPS constellation. The 
spacecraft is inertially pointing on an orbit whose average altitude is 575 km, its 
inclination is 57 degrees, and its eccentricity is 0.0015. The simulation generated 
from the RXTE ephemeris consists of 12 hours of sensor and reaction wheel data. 
After 6 hours, an attitude maneuver is simulated about the y axis at a rate of 0.2 
deg/sec for 6.63 minutes (giving a total maneuver of 80 degrees). Two GPS antennas 
are aligned along the spacecraft x and y body axes and one at the origin, pointing in 
the -z body axis direction. A magnetometer is aligned along the spacecraft body 
axes. The simulated magnetometer and GPS phase and range data contain added 
white noise, with standard deviations of 1 milliGauss, 1 degree, and 2 m, 
respectively. All the data are output at a rate of 2 Hz. The initial errors between the 
true state and the a priori filter state are given in Table I, except for the states b and z. 
The states b and z are initialized with the first magnetometer and GPS phase 
measurements, respectively. 

TABLE 1: INITIAL A PRIORI ERRORS IN THE FILTER STATE 


STATE 

Position 

Velocity 

Quaternion 

Rate 


ERROR 

llOOkm(RSS) 

0.6 km (RSS) 

1 1 1 deg rotation error 
0.001 deg/sec/axis 


RESULTS 

Figures 1 and 2 show the RSS of the angular velocity errors for measurements 
from one and two GPS satellites, respectively, during the first 150 minutes of the data 
span. In both cases the estimation approach is used (as applied to the GPS phase 
measurements, recalling that the estimation approach is applied to the 
magnetometer). Here it can be seen that the use of two GPS satellites shortens the 
convergence time. After the transients have settled there is little difference between 
the two cases. The RSS rate errors using the derivative approach are slightly noisier, 
but show similar behavior to the estimation approach for both the one GPS and two 
GPS satellite scenarios. The jumps in the RSS errors within the first 50 minutes 
occur when there is a change to a different GPS satellite. 



FIGURE 1 . RSS Angular Velocity Error During First 150 Minutes, Measurements 

from One GPS Satellite 



FIGURE 2. RSS Angular Velocity Error During First 150 Minutes, Measurements 

from Two GPS Satellites 

Figures 3 and 4 show the estimated y angular velocity for the estimation and the 
derivative approaches, respectively, each using measurements from one GPS 
satellite. The time span in the figures is during the attitude maneuver, which occurs 
after 6 hours. While in Figures 1 and 2 we compared the effect of one and two GPS 
satellites (because the difference between the estimation and derivative approaches 
was negligible), in Figures 3 and 4 we compare the two approaches because here the 
difference stems from the approaches and the difference between using one or two 
GPS satellites is negligible. In both cases, the filter follows the maneuver. The other 





two axes, namely the x and z axes, stay approximately at 0 deg/sec, similar to the 
estimated y rate both before and after the maneuver. 



time (min) 


FIGURE 3. Estimated Angular Velocity About Y Axis (Estimation Approach) 



time (min) 


FIGURE 4. Estimated Angular Velocity About Y Axis (Derivative Approach) 

Figures 5 and 6 show the RSS attitude error during the first 150 minutes, using 
measurements from one and two GPS satellites, respectively, using the estimation 
approach. Again, the use of two GPS satellites shortens the convergence time. After 
the filter has converged the results are similar with final average errors of 
approximately 0.3 degrees (computed over the last 30 minutes). The derivative 
approach results are similar to the estimation approach results. 





FIGURE 5. RSS Attitude Error During First 150 Minutes, Measurements from One 

GPS 



FIGURE 6. RSS Attitude Error During First 150 Minutes, Measurements from Two 

GPS 

Finally, Figures 7 and 8 show the RSS position errors for measurements from one 
and two GPS satellites, respectively, using the estimation approach. Here, the 
difference between using one or two GPS satellites is the most significant (the same 
is true for the velocity errors which have plots similar in shape to the position errors). 
After 12 hours, the one GPS case is still converging, with peak errors of about 4.5 km 
and an average RSS over the last 30 minutes of 0.8 km. The spikes in the errors 
occur due to a change in the GPS satellite and there is also an error from the magnetic 




field that occurs at the orbital period (approximately 90 minutes). In the case of 
using measurements from two GPS satellites, the effects of the changing GPS 
satellites are reduced. After 12 hours, the peak errors are reduced to less than 0.3 km 
with an average over the last 30 minutes of 30 m. The results for the final average 
RSS velocity errors are 6.6 cm/sec using two GPS versus 100 cm/sec when using one 
GPS. 



FIGURE 7. RSS Position Error, Measurements from One GPS 



FIGURE 8. RSS Position Error, Measurements from Two GPS 
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APPENDIX - Development of the GPS Line of Sight Derivative 

The position vector from the spacecraft to a GPS satellite is given as 

R s = Rsv-R (A.l) 

where R sv = the position vector of the GPS satellite from the center of Earth 
R = the position vector of the satellite from the center of Earth 

The GPS line of sight unit vector is then given as 



Taking the derivative of s results in 


s =■ 


dt 


v R *iy 


R s +r^ R s 

| R s| 


(A.3) 


The first term on the right hand side of (A.3) can be written as 


£ 

dt 


f 1 A 


VI R siy 


Rs-Vs 

|RJ 3 


(A. 4) 


where: V s = V sv - V = R s = the relative velocity of the GPS satellite with respect to 
the spacecraft 


Substituting (A.4) into (A. 3) results in 

^ R.-V. 


1 


s = — - 


R. 


3 Rs + |R s | Vs 


(A.5) 


The GPS position and velocity vectors R sv and V sv , respectively, are provided. 
However, only estimates of the spacecraft position and velocity vectors R and V, 
respectively, are available. Resolving all the vectors in inertial coordinates, and 
substituting in the estimated values results in the following expression for the 
derivative of the GPS line of sight vector implemented into the EKF algorithm. 


S,1 

-R,i+- 

— 

R„ 


R„ 


V S, 


(A. 6) 


where the A symbol refers to estimated values, and the ‘I’ refers to inertial 
coordinates. 


